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Abstract 

In this work, we investigate how and to which extent a quantum system can be driven along a 
prescribed path in Hilbert space by a suitably shaped laser pulse. To calculate the optimal, i.e., the 
variationally best pulse, a properly defined functional is maximized. This leads to a monotonically 
convergent algorithm which is computationally not more expensive than the standard optimal- 
control techniques to push a system, without specifying the path, from a given initial to a given 
final state. The method is successfully applied to drive the time-dependent density along a given 
trajectory in real space and to control the time-dependent occupation numbers of a two-level system 
and of a one-dimensional model for the hydrogen atom. 

PACS numbers: 42.50.Ct,32.80.Qk,02.60.Pn 
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I. INTRODUCTION 



Given a quantum-mechanical system, which laser pulse is able to drive the system from 
state A to state B in a finite time-interval? Which laser pulse maximizes the density in a 
certain given region in space by the end of the pulse? Questions of this kind are addressed 
by optimal control theory (OCT) in the context of nonrelativistic quantum mechanics. 

OCT as a field of mathematics dates back to the late 1950s and is widely applied in 
engineering. One of the most famous examples in engineering is the reentry problem of a 
space vehicle into the earth's atmospherefsee, e.g., Q]). The application of OCT to quantum 
mechanics started in the 1980s 3, [JO,^ ]. Due to the enormous progress in the shaping of 

became within reach. Experiments using 
closed loop learning (CLL) pj delivered highly convincing results 111. 

Calculated pulse shapes may be employed directly in the experimental setup, e.g., as an 
initial guess for CLL genetic algorithms. Perhaps more important, the theory can be used 
to decipher the control mechanism embedded in the experimental pulse shapes jl^ . 



The optimal control schemes 
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14 I15I employed so far in theoretical simulations and 
the experimental applications have been designed to reach a predefined target at the end of a 
finite time-interval. Little is known about controlling the path the quantum system takes to 
the desired target, i.e., controlling the trajectory in real space or in quantum number space. 
To our knowledge, three methods have been proposed so far: A fourth-order Euler-Lagrange 
equation to determine the envelope of the control field has been derived in Ref. [lfj]. This 
work, however, is restricted to very simple quantum systems. 

Another very elegant method, known as tracking, has been proposed by the authors of 
17] and Ref. Q|. Despite its tremendous success, this method bears an intrinsic 



Ref. 



difficulty: One has to prescribe a path that is controllable, otherwise singularities in the 
field appear, because of the one-to-one correspondence between the control field and the 
given trajectory. In practice, this may require a lot of intuition. 

The third method is an optimal control scheme for time-dependent targ ets jlfll ]. Basically, 
control schemes for time- independent targets [f^-bol] with an extension 



it combines optima 
to Liouville-space 



2l(. The schemes are generalized by introducing two new parameters 



like in Ref. |15] and then extended to include time-dependent targets. The new method is 
monotonically convergent and, in contrast to tracking, does not require a large amount of 
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intuition, i.e., choosing controllable obj ectives. Furthermore, the method is not restricted 
to two-level systems. While Ref. [19( presents the monotonically convergent algorithm, 
the full power of the method has not been exploited as yet. The challenge is the control 
of a truly time- dependent target represented by a positive-semidefinite, explicitly time- 
dependent operator. In Sec. II, we describe the general theory along with some examples of 
such operators. In particular, we discuss the control of occupation numbers in time and the 
indirect optimization of the dipole operator. The iterative procedure and some numerical 
details are explained in Sec. III. The results are presented in Sec. IV. 

II. THEORY 

We consider an electron in an external potential V(r) under the influence of a laser field. 
Given an initial state \I>(r, 0) = 0(r), the time evolution of the electron is described by the 
time-dependent Schrodinger equation with the laser field modeled in the dipole approxima- 
tion (length gauge), 

ij^M) = H*(r,t), (1) 
H = H -fie(t), (2) 
H = f + V (3) 

(atomic units are used throughout: h — m — e — 1). Here, fi = (£i x , £i y , £l z ) is the dipole 
operator and e(t) = (e x (t),e y (t),e z (t)) is the time- dependent electric field. The kinetic 
energy operator is T = — V 2 /2. 

Our goal is to control the time evolution of the electron by the external field in a way 
that the time- averaged expectation value of the target operator 0{t) is maximized. Mathe- 
matically, this goal corresponds to maximizing the functional 

Jim = I £dt (*(t)\d(t)Mt)), (4) 

where 0(t) is assumed to be positive-semidefinite. 

We want to keep the meaning of the operator 0(t) as general as possible at this point. 
A few examples will be discussed at the end of this section. 

Let us define 

3(t) = dx{t) + 2TS(t-T)d 2 , (5) 



so we can also include targets in our formulation that only depend on the final time T 

The functional J\\^\ will be maximized subject to a number of physical constraints. The 
idea is to cast also these constraints in a suitable functional form and then calculate the 
total variation. Subsequently, we set the total variation to zero and find a set of coupled 
partial differential equations P, Q. The solution of these equations will yield the desired 
laser field e(t). 

In more detail, optimizing J\ may possibly lead to fields with very high, or even infinite, 
total intensity. In order to avoid these strong fields, we include an additional term in the 
functional which penalizes the total energy of the field, 

J 2 [e] = -a [ dte 2 {t). (6) 
Jo 

The penalty factor a is a positive parameter used to weight this part of the functional against 
the other parts. 

The constraint that the electron's wave-function has to fulfill the time-dependent 
Schrodinger equation is expressed by 

J 3 [e,*, X ] = -23 J dt l^{t)\[id t -H)\^{t)) (7) 

with a Lagrange multiplier x(r, t). \l/(r, t) is the wave function driven by the laser field e(t). 
The Lagrange functional has the form 

J[ X , e] = Jif*] + J 2 [e] + J 3 [ X , e]. (8) 

Setting the variations of the functional with respect to \i ^ ' ■> an d e independently to zero 
yields 

ae .(t) = -QtxWfaMt)), J=x,y,z (9) 
= (id t -H^(r,t), (10) 
*(r,0) = 0(r), (11) 
(ja t -H)x{T,t) + ~d 1 {t)*{T,t) = 

i(x(r,t)-d 2 (tWr,t))6(t-T). (12) 
Equation Q determines the field from the wave function \l/(r, t) and the Lagrange multiplier 
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Equation ([10)1 is a time-dependent Schrodinger equation for \I/(r, t) starting from a given 
initial state 0(r) and driven by the field e(t). If we require the Lagrange multiplier x( r )^) 
to be continuous, we can solve the following two equations instead of Eq. (f*H?|) : 

(id t - H) X (r, t) = -^^(^(r, t), (13) 
X (r,T) = d 2 *(r,T), (14) 

To show this we integrate over Eq. (|12|) . 



lim 



jf + dti(x{r,t)-d 2 (t)*(r,t)y{t-T). (15) 



The left-hand side of Eq. f!15|) is because the integrand is a continuous function. It follows 
that also the right side must be 0, which implies Eq. (j!4)l . From Eqs. (114)1 and (fT2|) then 
follows Eq. PI) . 

Hence, the Lagrange multiplier satisfies an inhomogeneous Schrodinger equation with an 
initial condition at t = T. Its solution can be formally written as 

X (r,t) = Ul ]X (r,t ) - 1 /"dr ^(0!(r) *(r,r)), (16) 

where L^* is the time-evolution operator defined as C/* = Texp ^— i J^rft' H(t')^. 

The set of equations that we need to solve is now complete: Eqs. ©, (fTU|) . (|TT|) . (fl~3^) . 
and ()14|) . To find an optimal field e(t) from these equations we use an iterative algorithm 
which is discussed in the next section. 

In principle, we are not restricted to a single particle. The derivation and the algo- 
rithm can be generalized to many-particle systems, but except for a few model systems the 
numerical solution of the many-particle time-dependent Schrodinger equation is not feasible. 

We conclude this section with a few examples for the target operator 0(t). 

a. Final-time control. Since our approach is a generalization of the traditional optimal 



control formulation given in 
as a limiting case by setting 
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we first observe that the latter is trivially recovered 



Oi(t)=0, 2 = P=|$ / )($ / |. (17) 
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Here $/ represents the target final state which the propagated wave function \I/(r, t) is 
supposed to reach at time T. In this case, the target functional reduces to 

J x = (*{T)\P\*{T)) = \(^(T)\§f)\ 2 . (18) 



The target operator may also be local, as pointed out in Ref. [20|. If we choose 0\{t) = 
and 2 = 5{v — r ) (the density operator), we can maximize the probability density in r at 
t = T, 

J 1 = Jdx (y(T)\d 2 \y(T)) = n(T ,T). (19) 

Numerically, the 5 function can be approximated by a sharp Gaussian function. 

b. Wave-function-follower: The most ambitious goal is to find the pulse that forces 
the system to follow a predefined wave function $(r, t). If we choose 

o,{t) = mt))(mi (20) 

2 = 0, (21) 

the maximization of the time-averaged expectation value of Of with respect to the field 
e(t) becomes almost equivalent to the inversion of the Schrodinger equation, i.e., for a given 
function $(r, t) we find the field e(t) so that the propagated wave function \I/(r, t) comes 
as close as possible to the target $(r, t) in the space of admissible control fields. We can 
apply this method to the control of time-dependent occupation numbers, if we choose the 
time-dependent target to be 

\$(t)) = a (t)e- i£ot \0) + at^e-^ll) + a 2 (t)e- iS ^\2) + . . . , (22) 
H \n) = £ n \n), (23) 

o^t) = m))m)\. (24) 

The functions |a (t)| 2 , |ai(t)| 2 , |a 2 (t)| 2 , . . . are the predefined time-dependent level oc- 
cupations which the optimal laser pulse will try to achieve. In general, the functions 
ao(i), ai(t), a 2 (t), . . . can be complex, but as demonstrated in Sec. IV, real functions are suf- 
ficient in this case to control the occupations in time. For example, if in a two-level-system 
the occupation is supposed to oscillate with frequency Q we could choose a (t) = cos (fit) 
and, by normalization, a\(t) = sin(f2t). This defines the time-dependent target operator 
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c. Moving density. The operator used in Eq. (|19|) can be generalized to 



Oi(t) = 6(r-r (t)), (25) 



1 <- T 



J i = fj dt (*(*)i 5 ( r - r o(*))|*(*)) 

-T 



= i jT dtn(r (t),t). (26) 
Jx is maximal if the field is able to maximize the density along the trajectory r (t). 

III. ALGORITHM AND NUMERICAL DETAILS 

Equipped with the control equations 0, (|10p. and (fT3*j) we have to a find an algorithm to 
solve these equations for e(t). In the following, we describe such a scheme which is similar 
toRef. |l9j: 

zero-thstep: *W(0) ^> tfW(T) 
kthstep: [*W(T) ^ #W(0)] 

% «(T) ^ x «(0) 

[ X W(0) X (*)(T)] 



[*(*)( ) i!4 #«(T)] 
(0) ^ ^( fc+1 )(T). 



The laser fields used for the propagation are given by 

if\t) = (1 - r,)ef\t) - ^( X ^(tm¥ k \t)), (27) 

a 



e- ' '(f) = (l-T^W-^X^WlAii^H*)) J=*,V,*. (28) 



J J a 

The initial conditions in every iteration step are 

*(r,0) = 0(r), (29) 
X {r,T) = O a *(r,T). (30) 

The propagations in brackets are necessary only if one wants to avoid the storage of the time- 
dependent wave function and Lagrange multiplier. Note that the main difference between 
this iteration and the schemes used in is that one needs to know the time- dependent 
wave function ^f(r,t) to solve the inhomogeneous equation (JTSj) for the Lagrange multiplier 



x(r, t). Depending on the operator Oi(t), if the inhomogeneity is space- and time-dependent, 
this may require an additional time propagation. 

The choice of 77 and 7 completes the algorithm. 7 = 1 and 77 = 1 correspond to the 
algorithm suggested in |20j, while the choice 7 = 1 and rj = is analogous to the method 
used in j^J with a direct feedback of ^^(r, t). Further choices are discussed in Refs. [l^Q- 
In the following, we demonstrate the application of our algorithm to two different kinds of 
time-dependent targets, namely the control of occupation numbers and the control of a local 
operator. The first example is a two-level system consisting of states |0), |1) with a resonance 
frequency of 0% = u — uii = 0.395 and the di pol e matrix element P i — (1|A|0) = 1.05. 
The second system is a ID model for hydrogen 2^|, that has a "soft" Coulomb potential, 

V(x) = - T 4=T ( 31 ) 
Vr + 1 

This type of potential has been used extensively to gain qualitative insights in the behavior 



of atoms in strong laser pulses 



23, 



24 ]. 



The parameters Uqi and P01 of the two-level system are chosen to be identical with the 
lowest excitation energy and the corresponding dipole matrix element of ID hydrogen. 

The solution of the optimal control Eqs. Q, (fTUj) . (fT3*j) . and flTIj) requires the integration 
of the time-dependent Schrodinger equation with and without inhomogeneity. 
In the case of the two-level system, one may diagonalize the Hamilton operator analytically 
and therefore calculate the infinitesimal time-evolution operator directly. 

The time-dependent Schrodinger equation for the ID hydrogen model is solved on a 
grid, where the infinitesimal time-evolution operator is approximated by the second-order 
split-operator technique (SPO) 2s| . 

-t+At 



/ rt+At \ 

U l t +At = Texpf-z dt' H(t')j 



exp(-^ f At) exp(-i V(t) At) 

exp{~f At) + 0{At 3 ). (32) 



8 



For the inhomogeneous Schrodinger equation (j!4j). the infinitesimal time evolution of x( x > 



where we found the above, lowest-order approximation of the integral to be sufficient. 
Following the scheme described above, one needs five propagations per iteration step (if we 
want to avoid storage of the wave function). Within the 2nd order split-operator scheme 

n 

each time step requires four fast Fourier transforms (FFT) 26], because we have to know 
the wave function and the Lagrange multiplier in real space at every time step to be able 
to evaluate the field from Eq. (|9"|) . This sums up to 2 * 10 6 FFTs per 10 5 time steps and 
iteration. In comparison, optimal control methods for time-independent targets [15[ require 
four propagations. 

Since our hydrogen model can experience ionization, we employ absorbing boundaries to 
take care of boundary effects (otherwise we will find the outgoing wave incoming from the 
opposite boundary due to the periodic boundaries introduced by the Fourier transform). 
The real-space wave function is multiplied with a mask function that falls off like cos 1 - 1 / 8 ) 
at the boundary in every time step. 

IV. RESULTS 

1. Occupation number control 

First, we present the results for the two- level system. The time-dependent target wave 
function is chosen as |$(t)) = ao(t)e~ i£ °*|0) + ai(t)e~ l£lt |1), where the coefficients ao(t) and 
ai(t) are real and satisfy a§(i) + a?(t) = 1. 6 X = |$(t))($(t)|, 6 2 = 0. 

With the parameters At = 0.01, a = 0.05 and the initial guess field eo(t) = 10 -4 , the 
algorithm converges to a final value of J± = 0.9995 with a difference between two consecutive 
values of the functional 5J {n ' n+1) < 10" 8 . Fi gure^ shows the numerical results for the time 
evolution of the occupation numbers [Fig. |l(a)| and the optimized field [Fig. |l(b)[ . Figure 



is given by, 



X (x,t + At) 
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|l(c)| illustrates the monotonic convergence of the functional J\ + J2. The agreement between 
the calculated occupation and the V-shaped target function, as shown in Fig. |l(a)| is quite 
remarkable. To illustrate the quality of results associated with different values of J\, we 
have plotted the occupation curves corresponding to J x = 0.90, 0.95, 0.99 in Fig. |l(a)| 
Somewhat surprisingly, even if we reach J x = 0.95, there is still a sizable difference between 
the curves. Figure |l(b)| shows the envelopes extracted from the laser fields corresponding 
to Ji = 0.90, 0.95 as well as the optimal field corresponding to J\ = 0.9995. The resonance 
frequency of the system is found within a few steps. Then the algorithm improves the 
envelope. Figure |l(c)| shows the typical convergence behavior: a rapid improvement of the 
functional in the first few steps, implying that the difference between two consecutive fields 
is large (|46|) . and a slower convergence for the later steps, meaning that the differences 
between two steps get smaller [see Fig. |l(b)| . 

The same problem was solved for the ID hydrogen model on a grid. We found J\ = 0.97 
with 5J {n > n+1) < 10~ 5 . The parameters were At = 0.005, a = 1.5, and the initial choice 
for the field was again eo(t) = 10~ 4 . From our experience with the two-level system, we 
expect that the correspondence between the target curve and the optimized curve will not 
be perfect. 

The corresponding numerical results are shown in Fig. El Note, that the occupation of 
higher levels is negligible and that ionization is less than 0.2%. The field is in the weak 
response regime. 

n 

Tracking versus optimal control. Zhu and Rabitz showed [Tj| how the exact field neces- 
sary to follow a given trajectory can be determined by means of Ehrenfest's theorem. The 
exact field, however, may have singularities, i.e., the prescribed trajectory is not controllable 
with a smooth field. If, like in the next example, the target occupation curves bl(t),bl(t) 
consist of step functions, the exact field must have 5 peaks and, as a consequence, the 
tracking method cannot produce any useful results. The optimal control approach followed 
in this paper finds the best compromise between field energy and overlap with the target, 
yielding reasonable results such as those shown in Fig. EJ At times where becomes 
discontinuous, the field has intense pulses (see Fig |3(b)| ) consisting of only a few oscillations 
with the resonance frequency. 

The time-dependent occupation numbers in [Fig. |3(a)| deviate slightly from the target 
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FIG. 1: Target function and calculated occupation numbers for J\ = 0.90,0.95,0.99,0.9995 (a), 
optimal field J% = 0.9995 and extracted envelopes for J\ = 0.90,0.95 (b), and the value of the 
functional J\ + J2 (c) . 



curve. They are "washed out" at the discontinuity points of the target curve. For larger 
values of the penalty factor we notice that this broadening of the steps is even more pro- 
nounced (see Fig |4(a)| ). In this case, the width of the pulse envelope becomes broader and 
the maximum field strength lower [Fig. |4(b)| . 
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FIG. 2: Target curves ao(t),ai(t) and calculated occupation numbers (a), optimal field (b), and 
the value of J\ and J\ + J2 (c). 

2. Local operator 

A very important quantity to control is the time-dependent dipole moment. This quantity, 
however, cannot be accessed directly because the dipole operator is not positive semidefinite. 
As an alternative, we choose the time-dependent density operator, 

0(t) = 5(x-r(t)) (33) 
e-lr-rWf* . (34) 
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FIG. 3: Target curves bo(t),bi(t) and calculated occupation numbers (a), optimized field (b), and 
functional J\ + J2 (c) for a = 0.05. 

Intuitively, the dipole moment will roughly follow the curve described by r(t) since r(t) 
governs the movement of the density. 

In the actual computations, we approximate the 5 function by a sharp Gaussian (|34j). 

To test this idea, we first have to choose a reasonable function r(t). For this purpose we 
solve the time-dependent Schrodfnger equation for the ID hydrogen model with a given laser 
field e(i). From the resulting wave function, we calculate the time-dependent expectation 
value r(t) = (x)(t). With the function r(t) we then build the target operator (|3~4"j) and start 
our optimization with the initial guess eo(t) = 1CT 4 . 
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FIG. 4: The influence of the penalty factor a: occupation numbers (a) and optimal fields (b). 
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FIG. 5: (x) calculated with target and optimized wave function (a), optimal field (b) (after 1000 
iterations) for parameters a = 0.5, At = 0.005, eo = 10 -3 . 



Figure |5(a)| shows that the expectation value {x) opt calculated with the optimal field 
e opt (t) follows the target r(t) rather closely. We do not obtain a perfect correspondence 
between (x)(t) and r(t), but the results clearly demonstrate that the algorithm also works 
for this type of target and, hence, that the indirect approach to control the dipole moment 
is appropriate. 

As proven in the Appendix, it is also possible to optimize functionals of the type J\ = 
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J Q T dt(^f (t)\0(t)\^f (t)) , n > 1. Since our integrand l\ is < 1, the effect will be that J\ carries 
h 

less weight in the optimization, i.e., the algorithm will try to decrease the field energy. Hence 
we expect the similarity between the target trajectory r(t) and the calculated expectation 
value (x) opt to be less for increasing n. The results for n = 2, 3, 4 are shown in Fig. |H1 If we 
build a target functional with the integrand I\ > 1 we will find the opposite effect, i.e., J\ 
will become more important than before. This demonstrates that the parameter n provides 
a new handle (in addition to the penalty factor a) to shift the relative importance of J\ 
versus J 2- 

0.2 
0.1 

3 Q, 

A 
X 
V 

-0.1 
-0.2 

50 100 150 200 250 300 
time (a.u.) 

FIG. 6: Comparison for the expectation value (x) with the target trajectory (squares) for different 
exponents n. 

V. CONCLUSION 

This work deals with the quantum control of time-dependent targets. In Sec. II, we 
presented explicit examples of positive-semidefinite target operators designed for the control 
of time-dependent occupations and of the time-dependent dipole moment. We then applied 
these operators to control the time evolution of a simple two-level system and of a grid model 
(ID hydrogen). In each case, we find a continuous increase of the value of the functional 
J\-\- The improvement in the first iteration steps is quite strong, while it takes a large 
number of iterations to converge the last few percentages. The results also show that a large 
number of iterations is required to reach perfect agreement with the target trajectories. In 
the Appendix we prove that by exponentiating the integrand with a positive integer n > 1, 
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the iteration still converges monotonically. The functional constructed in this way contains 
two parameters, a and n, that allow one to fine-tune the relative importance of J\ and J 2 . 
To summarize, we have demonstrated in this work that the quantum control of genuinely 
time-dependent targets is feasible. In particular, the successful control of the dipole moment 
in time may open new avenues to optimize high-harmonic generation, which is extremely 
important for shaping attosecond laser pulses. Work along these lines is in progress. 
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We want to show that the same iteration will converge monotonically also for a functional 
of the type 

h = [ T dt(*(t)\d(t)\y(t)y, 

Jo 

where n > 1, n G N. The equation for the Lagrange multiplier then has the form 

(id t - H) X (r, t) = -ni^^ld^m^dm^t), (35) 
X (r,T) = 0. (36) 

Now consider a and 6 as real positive numbers. Then 

a n_ b n = nh n-^ a _ fc) + + _ _ n& «-l Q . (37) 

N v ' 

=A 

Next we show that A = a n + (n — 1)6™ — nb n ~ l a is never negative. Defining a = b + 5, 
5 G [—6, oo), we distinguish between two cases. 
Case I: 5 G [0, oo), 

a n + (n - 1)6™ - nb n - l a 
= (b + 5) n + (n- 1)6™ - nb n -\b + 5) 
= 6™ + nb^'S + ■■■ + nS^b 
+5 n + (n- 1)6™ - nb^S 

= j-—^b n - 2 5 2 + --- + n5 n ~ 1 b + 5 n >0. (38) 
16 



Case II: 5 G [-6,0), 



{b + 5) n +{n- l)b n - nb n -\b + 5) 
5\ n 



>o 



1 + 



+n — 1 — n — n- . (39) 
o 

To evaluate Case II, we introduce x = S/b, x G [—1, 0) with, 

f(x) = (l + x) n -nx -1, 
y := x + 1, 

y e [o,i), 

f(y) = y n -ny + n-l, 
/(0) = n-l, 

/(I) = o, 

f'(y) = n^-l) 

<o 

/(y)G[0,n-l] 

=► /(v) > o. 

Since /'(y) < 0, the function / must decrease monotonically from n — 1 > to 0, so it 
cannot become negative. In conclusion, 



a n -b n = nb n -\a -b)+a n +(n- l)b n - n6 n_1 a . 



(40) 



A>0 



The deviation in J between two consecutive steps is, 



§j(k+i,k) = j(fc+i) _ j(k) 

= J T dt (^¥ k+1 \t)\d{t)\¥ k+1 \t)) n 
- {¥ k \t)\d(t)\¥ k \t)) n 

+ a[e^(t)] 2 -a[^ k+1 \t)] 2 ). 
If we identify a(t) = (^ k+1 \t)\d(t)\^ k+1 \t)) and b(t) = (& k \t)\d(t)\^ k \t)) and use Eq. 
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6J^ k+1 ^ = £dt (^A(t) + a[^ k \t)} 2 - a [ e ( fc +D(t)] 2 

+ n(¥ k \t)\d(t)\¥ k \t)) n - 1 

\¥ k+1 \t)\d(t)\¥ k+1) (t)) - (¥ k) (t)\d(t)\¥ k) (t)) 



where we have separated the positive term A(t) = a n (t) + (n — l)b n (t) — nb n 1 (t)a(t). 

We define B(t) = n(*W(*)|0(t)|*W(*)) n - 1 (<y*( fc+1,fe )(t)|0(t)|<y*f fc+1 ' fc >(t)) > and rewrite 
SJ {k+i,k) ag 



5.j( k+1 >V = Jdt (A(t) + B(t) + a [e {k \t)Y - a [e (fc+1) (*)] 

+ (t) \3(t) |^ (fc) (t)) n - 1 23?(^( fc ) (t) |6(t) < fe+1 - fe >(t)) 

where S^ {k+1 ' k \r,t) = & k+1 \r,t) - #W(r,t). We use the equation for the Lagrange mul- 
tiplier (}35j) and obtain 

5J (fe+1 ' fc) = jdt (A(t)+B(t) 

+ a[^ k \t)] 2 - a [^ k+1 \t)} 2 

+ 2&(- (d t + iH^ k Ax {k \t)\5¥ k+1 ' k \t)^, (41) 



where = H — fie^ k \ 



= Jdt (A(t) + B(t) 

+ a [ e ( k \t)} 2 - a [e^\t)} 2 

+ 2$(x {k) (t)\ (id t -H^ \6¥ k+1 ' k \t) 

+ 2K( X («)l^ (Jtfl '* ) (*)>lo- ( 42 ) 







For the last term in Eq. (g2D we used the fact that S^ k+1 ' k \r, 0) = since the initial state 
for the wave function is fixed and x(r, T) = because of Eq. (j3SJ). 

We use the time- dependent Schrodinger equation (jlUj) for and \fr( k+1 \ where = 
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#o-£e (fc) , 

(id t -H {k) ^j 5^ {k+hk \r,t) 
= (id t -H {k) ^j (*<* +1 >(r,t)-*W(r,f)) 
= (#( fc+1 )- J ffW)^( fc+1 )(r,t) 

_(#(*)_#(*>) ¥ k \r,t) 
= _ (e^ k+1 \t) -? (fc) (t)) A(r)^ (fc+1) M) 

+ (>)( t ) _e< fc) (t)) £( r )^ fc )(r,t). (43) 

Consequently, the change in the Lagrange functional becomes 

5J (fc+1 ' fe) = £dt (A{t) + B{t) + {5¥ k+1 ' k \t)\d{t)\8¥ k+1 ' k \t)) 

- a ([e^)] 2 -[ e «(t)] 2 ) 

-23(x (fc) (t)\fi\¥ k+1) (*)) (e (fc+1) (t) - ? (fe) (t)) 

+23( X (/c) WlAl^ (fc) W)(e (/£) W -^(i)))- (44) 
Finally, using Eqs. {2Zj) and we find 

eft (A(t) + B(t) + (^ (fe+1 ' fe) (t)|6(t)|^ (fc+1>fe) (t)) 



o \ 

I 2 . r mi2 



+ 2^ +1 )(t) - (1 - 7 )e (fc) (t))^(e (fc+1) (t) -^(t)^ 
- 2^(t)-(l-i 7 )eW(t))^e«(t)-i<«(t))) (45) 

dt \A(t) + B(t) + (5¥ k+1 ' k \t)\d(t)\5¥ k+1 ' k \t)) 



o 

+ a(^-l) (e( fc+1 )(f)-6«(t)) ; 



+ «Q-l)(6«(*)-^)) a Y (46) 

For ?7, 7 G [0, 2] (similar to Q and Q), the iteration converges monotonically, i.e., 
§j(k+i,k) y q This iteration converges monotonically and quadratically in terms of the 
field deviations between two iterations. 
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We emphasize that this proof is true only if the solution of the time-dependent 
(in)homogeneous Schrodinger equation is exact in each step. Numerical implementations 
are of course always approximate and, as a consequence, it may happen that the value of 
the functional J decreases. This, on the other hand provides a test of the accuracy of the 
propagation method. 
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